#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#define YASP
// #define DEBUG
#include "includes/need-a-name.hh"
using namespace Dune;
int main(int argc, char** argv)
{
  try{
    // initialize MPI
    MPIHelper& helper = MPIHelper::instance(argc, argv);
    std::cout << "Need-A-Name, a dg mhd code for multi-scale multi-physical problems." << std::endl;
#ifdef DEBUG    
    if(MPIHelper::isFake)
      std::cout<< "This is a sequential program." << std::endl;
    else
      std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
        <<" processes!"<<std::endl;
#endif
    // create a project
    // ProjectBase* aproject = create_project(constantAdvection_1D_YASP);
    // ProjectBase* aproject = create_project(soundWave_1D_YASP); 
    // ProjectBase* aproject = create_project(euler_1D_YASP);
    // ProjectBase* aproject = create_project(mhd_1D_YASP);
    // ProjectBase* aproject = create_project(fv_2D_YASP);
    ProjectBase* aproject = create_project(fv_3D_YASP);
    // ProjectBase* aproject = create_project(dg_1D_YASP);
    aproject->set_project_signature();
    std::cout<< "Project signature: " << aproject->project_signature << std::endl;
    aproject->load_subprojects();

    // subproject loop
    for(int subproj=0; subproj<aproject->get_subproject_number(); subproj++){
#ifdef DEBUG
      std::cout <<"subproject index; "<< subproj << std::endl;
#endif
      SubprojectBase* asubproject = aproject->get_subproject(subproj);
      asubproject->set_subproject_signature();
      std::cout<< "Subproject signature: " << asubproject->subproject_signature << std::endl;
      // prepare memory
      memory::variable.erase(memory::variable.begin(),memory::variable.end());
      std::vector<std::vector<double> >(memory::variable).swap(memory::variable);
      memory::variable_tmp.erase(memory::variable_tmp.begin(),memory::variable_tmp.end());
      std::vector<std::vector<double> >(memory::variable_tmp).swap(memory::variable_tmp);      
      memory::variablename.erase(memory::variablename.begin(),memory::variablename.end());
      std::vector<std::string>(memory::variablename).swap(memory::variablename);
      memory::limit_slope.erase(memory::limit_slope.begin(),memory::limit_slope.end());
      std::vector<std::vector<std::vector<double> > >(memory::limit_slope).swap(memory::limit_slope);
      memory::update.erase(memory::update.begin(),memory::update.end());
      // std::vector<std::vector<double> >(memory::update).swap(memory::update);
      std::vector<std::vector<std::vector<double> > >(memory::update).swap(memory::update);
      memory::local_vmax.erase(memory::local_vmax.begin(),memory::local_vmax.end());
      std::vector<double>(memory::local_vmax).swap(memory::local_vmax);
      dataanalysis::auxvariablename.erase(dataanalysis::auxvariablename.begin(),dataanalysis::auxvariablename.end());
      std::vector<std::string>(dataanalysis::auxvariablename).swap(dataanalysis::auxvariablename);  
      // assign equation, initial-condition and boundary-driver
      asubproject->assign_equation();
      asubproject->assign_initialcondition();
      asubproject->assign_boundarycondition();      
      // set parameters
      aproject->set_common_parameters(subproj);
      asubproject->set_unique_parameters();
      asubproject->set_variable_name();
  
      // create grid
      // need to improve in later version
#ifdef YASP
      switch (setting::spacegrid){
       case setting::yasp:
	std::cout << "space grid: yasp" << std::endl;
	break;
	case setting::ug:
	std::cout << "space grid: ug" << std::endl;
	break;
      }
      std::cout<<"dimension: "<<gridsetting::dimension << std::endl;
      int overlapSize = 1;      
      using Grid = YaspGrid<gridsetting::dimension,EquidistantOffsetCoordinates<double,gridsetting::dimension>>;
      std::bitset<gridsetting::dimension> periodic("000");
      FieldVector<double,gridsetting::dimension> lowervalue;
      FieldVector<double,gridsetting::dimension> uppervalue;
      std::array<int,gridsetting::dimension> elementvalue;
      for(int tmpi=0;tmpi<gridsetting::dimension;tmpi++){
	lowervalue[tmpi]=gridsetting::lower[tmpi];
	uppervalue[tmpi]=gridsetting::upper[tmpi];
	elementvalue[tmpi]=gridsetting::element[tmpi];	
      }
      std::cout << "create calculation grid" << std::endl;            
      Grid* agrid= new Grid(lowervalue,uppervalue,elementvalue,periodic,overlapSize,MPI_COMM_WORLD);
#endif
#ifdef UG
      std::cout << "using UGGrid" << std::endl;
#endif 

      // gridview and attach memory
      std::cout << "attach memory" << std::endl;      
      using GridView = Grid::LeafGridView;
      GridView gridView = agrid->leafGridView();
      MultipleCodimMultipleGeomTypeMapper<GridView> mapper(gridView, mcmgElementLayout());
      memory::variable = std::vector<std::vector<double> >((memory::dgorder+1)*memory::phynumber+dataanalysis::auxnumber, std::vector<double>(mapper.size()));
      memory::variable_tmp = std::vector<std::vector<double> >(std::max(1,memory::dgorder)*memory::phynumber+dataanalysis::auxnumber, std::vector<double>(mapper.size()));      
      switch(setting::timeadvance){
      case setting::euler:
	memory::updatenumber=1;	
	break;
      case setting::rk2:
	memory::updatenumber=2;		
	break;	
      case setting::rk3:
	memory::updatenumber=3;	
	break;
      case setting::rk4:
	memory::updatenumber=4;
	break;
      }
      memory::update = std::vector<std::vector<std::vector<double> > >(memory::updatenumber,std::vector<std::vector<double> >(std::max(1,memory::dgorder)*memory::phynumber+dataanalysis::auxnumber,std::vector<double>(mapper.size())));
      memory::limit_slope = std::vector<std::vector<std::vector<double> > >(memory::phynumber,std::vector<std::vector<double> >(gridsetting::dimension,std::vector<double>(mapper.size())));      
      memory::local_vmax = std::vector<double>(mapper.size());
      
      // initialization, vtk writer, init-information  
      GoThroughEntity::go("initialization",asubproject,gridView,mapper);
      auto vtkWriter = std::make_shared<Dune::VTKWriter<GridView>>(gridView);
      VTKSequenceWriter<GridView> vtkSequenceWriter(vtkWriter, aproject->project_signature+"___"+asubproject->subproject_signature);
      for(int index=0; index<((memory::dgorder+1)*memory::phynumber+dataanalysis::auxnumber); index++)
	vtkWriter->addCellData(memory::variable[index],memory::variablename[index]);
      std::cout << "output initial snapshot" << std::endl;            
      vtkSequenceWriter.write(timerule::current);
      
      // time loop
#ifdef DEBUG 
      std::cout << "timerule::current: " << timerule::current << " timerule::end: " << timerule::end << std::endl;
      // only test init
      timerule::end=0.0;
#endif
      std::cout << "time loop begin" << std::endl;
      while (timerule::current<=timerule::end && timerule::counter<=200){
	// std::cout << "into the time-loop"  << std::endl;
 	GoThroughEntity::go("dt",asubproject,gridView,mapper);
	GoThroughEntity::go("slope",asubproject,gridView,mapper);	
        // iteration information
	if(memory::dgflag==0)
	  GoThroughEntity::go("fvevolve",asubproject,gridView,mapper);
	if(memory::dgflag==1){
	  GoThroughEntity::go("dgevolve_finiteElement",asubproject,gridView,mapper);
	  GoThroughEntity::go("dgevolve_calculateFlux",asubproject,gridView,mapper);
	  GoThroughEntity::go("dgevolve_combineWeight",asubproject,gridView,mapper);	  
	}
	timerule::current += timerule::dt;
        if(helper.rank()==0) {
          std::cout << "counter=" << timerule::counter << " time=" << timerule::current << " dt=" << timerule::dt << std::endl;
        }
        timerule::counter++;
        // output rule
	if (timerule::counter%5==0){
	  asubproject->dataanalysis();
          vtkSequenceWriter.write(timerule::current);
	  }
      }
      // relase memory
      memory::update.erase(memory::update.begin(),memory::update.end());
      std::vector<std::vector<std::vector<double> > >(memory::update).swap(memory::update);
      memory::limit_slope.erase(memory::limit_slope.begin(),memory::limit_slope.end());
      std::vector<std::vector<std::vector<double> > >(memory::limit_slope).swap(memory::limit_slope);
      memory::local_vmax.erase(memory::local_vmax.begin(),memory::local_vmax.end());
      std::vector<double>(memory::local_vmax).swap(memory::local_vmax);
      // relase grid
      delete agrid;
      agrid=NULL;
      std::cout<< "release grid"<<std::endl;
    }
  }
  catch (Exception &e){
    std::cerr << "Dune reported error: " << e << std::endl;
  }
  catch (...){
    std::cerr << "Unknown exception thrown!" << std::endl;
  }  
  return 0;  
}
